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Numerical simulations of rotating axisymmetric sunspots 
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ABSTRACT 

A numerical model of axisymmetric convection in the presence of a vertical magnetic 
flux bundle and rotation about the axis is presented. The model contains a compress- 
ible plasma described by the nonlinear MHD equations, with density and temperature 
gradients simulating the upper layer of the sun's convection zone. The solutions ex- 
hibit a central magnetic flux tube in a cylindrical numerical domain, with convection 
cells forming collar flows around the tube. When the numerical domain is rotated 
with a constant angular velocity, the plasma forms a Rankine vortex, with the plasma 
rotating as a rigid body where the magnetic field is strong, as in the flux tube, while 
experiencing sheared azimuthal flow in the surrounding convection cells, forming a 
free vortex. As a result, the azimuthal velocity component has its maximum value 
close to the outer edge of the flux tube. The azimuthal flow inside the magnetic flux 
tube and the vortex flow are prograde relative to the rotating cylindrical reference 
frame. A retrograde flow appears at the outer wall. The most significant convection 
cell outside the flux tube is the location for the maximum value of the azimuthal mag- 
netic field component. The azimuthal flow and magnetic structure are not generated 
spontaneously, but decay exponentially in the absence of any imposed rotation of the 
cylindrical domain. 
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1 INTRODUCTION 

Observations of sunspots rotating around their own axis, 
which is perpendicular to the plane of the photo spher e , have 
a long history throughout the twentieth century. faald (|l908l l 
and Evershed (1910) noticed the rotation as well as a vor- 
tex forming around the rotating sunspot. Observations in 
the photosphere and corona continued through the century, 
culmin ating in th e high resolution measurements of today. 
fSee iBrown et all (|2003") and references therein.) An excit- 
ing development during the last decade has been the abil- 
ity to measure the associated flow beneath the photosphere 
(|Gizon Birchl[2005l : IZhao Kosovichevll2003l l. 

There is a distinct radial profile as sociated with t he az - 
imuthal velocity of rotating sunspots. 'Bro wn et al.l (|2003l ) 
found that the umbra (in which the rotation axis resides) 
has small average azimuthal velocities, while the fastest ro- 
tation occurs at some point along the radial length of the 
penumbra. The rotation then tails away to a negligible value 
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outside of the sunspot. The peak azimuthal velocity in the 
penumbra can be m ore than double that inside the umbra. 
iBrown et al.l (|2003l l found suggestions of rotation outside 
some sunspots, but these observations ar e hampered by an 
ambiguous penumbral edge. In contrast, lYan Qui (|2007f ) 
observed a rotating sunspot where the maximum azimuthal 
velocity occurred inside the umbra. The rotation persisted in 
the penumbra and the area near the penumbra, with the an- 
gular velocity reducing as one moves radially away from the 
umbra. The surrounding area far removed from the penum- 
bra experienced a slow rotation in the opposite direction 
from that of the rotating sunspot. 

It is not clear if the di r ectio n of rotation has a hemi- 
spheric preference. iKnoskal (| 19751 ). and references therein, 
found that the majority of rotating sunspots in both hemi- 
spheres turn anticlockwise. However, iDing, Hong Wan j 
(|l987l l found a preference, with clockwise (anticlockwise) 
rotation predominantly in the southern (northern) hemi- 
sphere. This would suggest that the Sun's differential ro- 
tation associated with the global flow field has an influ- 
ence. The Coriolis force due to the Evershed flow field 
would cause rotation in the opp osite direction. How- 
ever, helioseismic observations (Giz on BirchI 20051)^ sup- 
ported by numerical results (|Hurlburt &: Rucklidgd l2000l : 
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iBotha. Rucklidge Hurlburtll2006h . show a converging hor- 
izontal flow below the Evershed flow, which leads to cyclonic 
vorticity and hence a possible contribution from the Coriolis 
force. 

The sni a ll sam ple of rotating sunspots studied by 
iBrown et alj (|2003l l suggests that younger sunspots rotate 
faster than older ones. However, it was difficult to judge the 
ages of the sunspots in the sample. The rotation rates are 
time dependent, with all rotation eventually decreasing with 
time. The peak rotation in the penumb ra can be anything u p 
to 3 degrees per hour, as observed bv iBrown et al.l (|2003l ). 
The same time evolut ion was observed in a rotating pore 
([Porotovic et al.l *2Q02) . This behaviour suggest that some of 
the rotation are caused by local events that are transitory in 
nature. This conclusion is strengthened by an observation of 
damped oscillatory motion, which had a maximum rotation 
of 3.5 degrees per hour (Kucera 1982). 

A possible mechanis m causing rotating sunspots is the 
rise of twisted flux ropes (jGibson et al.ll2004h . In this model 
the rotation of two flux rope poles is observed after the 
central horizontal portion of the flux rope has emerged 
through the photosphere. This implies the existence of two 
co-evolving sunspots of opposite magnetic polarity. This is 
generally not in evidence in the observations, mostly be- 
cause leading sunspots are often followed by a more diffuse 
opposite polarity. 

It was shown by iGopasvuk Gopasvukl (|2005l l that 
when the averaged velocity and magnetic fleld compo- 
nents are subtracted from observed sunspots, the result 
fits a lightly damped sinusoidal wave. This implies that as 
well as the motion described so far, sunspots also experi- 
ence torsional oscillat ions. Using the thin flux tube model, 
iMusielak et al.l (|2007l ) found that in a compressible, isother- 
mal fleld-free medium, linear torsional Alfven waves along 
the magnetic tube do not have a cutoff frequency. 

The rotation of sunspots has been linked to the 
formation of soft X-ray si g moid s as well as the 
erupt i on of flares (lAlexanderl l2006l : iRegnier &: Canfleldl 
'2006; "Tian k Alexa nderl 12006^ Numerical simulations by 
Gerrard et al. (2003) show that by adding a horizontal pho- 
tospheric flow to the rotation, the generation of flares is 
enhanced as both rotation and flow increase the complexity 
of the magnetic fleld. This is supported by the observation 
that flare activity is correlated to magnetic flux and kinetic 
vorticity (Mason et al. 2006). 

Evidence from helioseismic measurements show that the 
rotation of sunspots, as observed in the photosphere, extends 
into the deeper layers of the sun. Up to a depth of approx- 
imately 7 Mm vortical flow in the same direction as the 
rotation of the sunspot exists, while below 7 Mm a vortical 
flow opposite to the sunspot rotation direction is observed 
(lKosovichevll2002l : IZhao &: Kosovichevlliooi : iGizon BirchI 
l2005l ). However, it should be noted that helioseismic mea- 
surements are difficult and not always consistent. For ex- 
ample, up to a depth of 3 Mm a converging inflow is found 
when using p modes, whil e f mode measureme nts find only 
outflows down to 10 Mm (|Gizon k Birchll2005h . 

In this paper an axisymmetric model is used to sim- 
ulate rotation around a central magnetic flux bundle. The 
values of the physical parameters of the model are chosen to 
describe the solar convection zone from a depth of approxi- 
mately 500 km below the visible surface of the Sun to a depth 



of ap proximately 6000 km (|Botha, Rucklidge k HurlburtI 
I2OO6I I . The numerical domain is a cylinder with an aspect ra- 
tio of r = 3, i.e. one unit deep and a radial distance of three 
units. This implies that we are simulating magnetoconvec- 
tion on the supergranular scale. To generate azimuthal flow 
and magnetic fleld, the whole domain is rotated at a constant 
angular velocity. Strictly speaking this is not equivalent to 
simulating a pore or sunspot where only the magnetic flux 
bundle rotates. However, in spite of driving the azimuthal 
flow throughout the numerical domain, we f ind that the 
solution tends to conform to observations of iBrown et al.l 
f2003), where a maximum azimuthal flow occurs close to 
the magnetic flux bundle. This means that a vortical flow 
forms around the flux bundle while the plasma inside the 
bundle rotates as a solid body. This type of flow is formally 
described as a Rankine vortex. 

Our nu merical results may be compared to results found 
bv I Jones Gallowavl (| 19931 1, who studied a Boussinesq fluid 
in an axisymmetric cylinder. They imposed two types of 
boundary conditions: a stress-free outer wall as well as an 
external flow, implemented by rotating the outer wall of the 
flxed cylinder. A flux bundle formed at the central axis, with 
the maximum angular velocity occurring near the axis. For 
high magnetic field strength, described by the dimensionless 
Chandrasekhar number (Q), the fiux bundle broadened with 
the stress-free boundary conditions producing a reverse in 
azimuthal magnetic fiux near the axis, while the imposed 
external fiow produced a reverse azimuthal fiow near the 
axis. The reversal in the azimuthal magnetic field is ascribed 
to the conservation of angular momentum under stress-free 
conditions, while the fiow reversal obtained with the im- 
posed external fiow is ascribed to the working of the Lorentz 
force. 

In Section [2] the mathematical model and its numerical 
implementation is described. This is followed by the numer- 
ical results. When no rotation of the numerical domain is 
present (Section 13. 1|) . no azimuthal velocity and magnetic 
field are generated spontaneously: both quantities are small 
and decay exponentially with time. This case is useful to 
compare the rotating solutions against. Driven by the rotat- 
ing numerical domain, the solution settles into a time inde- 
pendent solution that shows rigid rotation of the plasma in 
the magnetic fiux tube and vortical rotation around it (Sec- 
tions [321 and [3]3]). This solution is robust and essentially 
stays the same when the magnetic field strength is increased 
(Section 13. 4p . the Prandtl number is lowered (Section 13. 5|) , 
as well as when the stratification is increased (Section 13. 7|) . 
The latter part of the numerical investigation explores the 
infiuence of the numerical domain on the solution. This we 
do by changing the bottom and outside boundary condi- 
tions (Sections 13.61 and 13. 8|) . We conclude the paper with a 
summary of the results. 



2 MODEL 

Partial differential equations describing compressible mag- 
netoconvection are solved in an axisymmetric cylindrical ge- 
ometry, using a numerical code developed for this purpose. 
A detaile d description of the tw o dime nsional (2D) model is 
given bv lHurlburt Rucklidgj (|2000h . Here we extend the 
model to 2.5D by including azimuthal components in addi- 
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tion to the radial and axial components. A constant angular 
velocity is added that introduces the Coriolis and centrifugal 
forces into the Navier- Stokes equation. 



2.1 Mathematical model 

The initial temperature and density profiles in the vertical 
(z) direction are given by 



T 
P 



(1) 

(2) 



The temperature and density are scaled so that they are 
equal to 1 at the top of the static atmosphere. The initial 
temperature gradient is given by ^, while m is the poly- 
tropic index. The equations for fully compressible, nonlinear 
axisymmetric magnetoconvection that we use are 
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The cylindrical reference frame is rotated about its axis at a 
constant angular velocity of = (d(/)/dt)z, which is respon- 
sible for the Coriolis and centrifugal terms in the Navier- 
Stokes equation JJ]). The vector potential A^f, gives the r 
and z components of the magnetic field while the azimuthal 
component is included explicitly, so that the magnetic field 
is given by 



B = V X {cj)A^) + (l)B^. 



(8) 



The velocity consists of three components, namely u = 
{ur,U(f),Uz), where refers to the azimuthal velocity rel- 
ative to the rotating reference frame. We also use the auxil- 
iary equations 



V-B = 0, 



j = VxB, 



(9) 



and the following notation: 7 is the ratio of specific heats; 
a the Prandtl number; ("o the magnetic diffusivity ratio at 
z = 0; and the Chandrasekhar number given by 
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(10) 



where d is the depth of the domain, ji the magnetic perme- 
ability, Tf the magnetic diffusivity and v the kinetic viscosity. 
The rate of strain tensor is given by 
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while the dimensionless thermal conductivity K is related 
to the Rayleigh number R in the following way: 



i? = r(m + l) 



(m+l)(7-l) 

7 



(1 + ^/2)' 



(12) 



R is du measure of the importance of buoyancy forces com- 
pared to viscous forces in the middle of the layer, and is 
used to drive the convection in the model. The addition of 
rotation adds an additional scaling which relates the convec- 
tive timescale to the rotational timescale. Following iGilmanI 
(1977) we express this ratio using the convective Rossby 
number defined as 



(13) 



iBrummell. Hurlburt Toomrd (|l996[ ) found that this pa- 
rameter, which can be evaluated from the control parame- 
ters, is typically close to that based on the traditional def- 
inition of the Rossby number, namely the ratio of the rms 
vorticity of the flow to the vorticity 2Q associated with the 
rotating cylinder. The effects of rotation are significant for 
Ro ~ 1, and dominate for Ro <^ 1. For supergranules and 
sunspot moats, Ro ^ 30. 

All the other symbols have their usual meaning. The 
physical quantities are dimensionless, with the length scaled 
proportional to the depth of the numerical domain, velocity 
scaled proportional to the sound speed at the top of the do- 
main and temperature, magnetic field, density, and pressure 
all scaled proportional to their initial values at the top of the 
numerical domain. These top initial values are radially uni- 
form and do not change throughout this paper, as discussed 
m Section [221 



2.2 Numerical implementation 

The computational domain is an axisymmetric cylinder of 
radius F, situated in the (r, z) plane so that 

O^r^F, O^z^l, (14) 

with ^ = at the top of the box ([Hurlburt Rucklidgd 
I2OOOI : iBotha. Rucklidge Hurlburtll2006l ). We require that 
all variables be sufficiently well-behaved at the axis (r = 0) 
and that the differential operators in the PDEs are non- 
singular. This implies that 



dp duz dT ^ 

Ur = = — — = — - = 

or or or 

A^^Br^B^^^^j^^ 0. (15) 

Terms like Ur/r, U(f)/r and P^/r are evaluated using 
I'Hopital's rule, while terms containing Urjr^ cancel alge- 
braically. 

The outside wall (r = F) is a slippery, impenetrable 
wall with no lateral heat flux across it (i.e. an insulator): 
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0, 



dz 



(18) 



g^=^r = ^=Br = B,=J,=0, -^ = -. (16) 

The magnetic potential has the value — V/2 at the out- 
side wall, which was chosen so that the initial vertical uni- 
form field satisfies Bz = 1. 

At the bottom boundary the magnetic field is vertical. 
The temperature T is chosen to be constant with value + 1 
from equation ([T]). The bottom boundary is impenetrable 
and stress free, i.e. 

„ „ dBz dur du^ 

^^ = ^^ = ^ = ^ = ^ = -^ = 0. (17) 

The top of the box is treated as impenetrable for the 
plasma, with a radiative temperature boundary condition 
given by Stefan's law: 

dUr _ dU(j) _ 

dz dz 

The in Stefan's law is the same as in ([T]), so that the 
equilibrium profile used as initial condition is not destroyed. 
The Stefan- Bolt zmann constant will enter (|18|) only when 
we dimensionalize it. The magnetic field is matched to a 
potential field on top of the numerical domain, dA(f,/dz = 
Mpot{A(j)), where Mpot is a linear operator, so that Br and 
Bz are continuous across the boundary. The potential field is 
solved by assuming an infinitely tall conducting cylinder of 
radius F above the domain, with the magnetic field becoming 
uniform as z ^ oo. A more d etailed description is given by 
iHurlburt k Ruckhdgd (|2000l ). No currents exist inside the 
potential field and consequently we choose jz = along the 
top of the box. From (|15p it then follows that i?^ = along 
the top boundary. 

One consequence of these boundary conditions is that 
no current escapes from the numerical domain: jr = on 
r = F and jz = on z = and 1. It follows that the 
boundaries do not provide any net vertical torques, and so 
do not contribute to changes in the vertical component of 
the total angular momentum, Lz. Nonetheless, Lz is not 
conserved in compressible convection: the Coriolis term can 
lead to changes in L^, for example, when mass is transported 
to larger distances from the axis. (Note that this does not 
occur in incompressible convection.) A consequence of this 
is that as the solution evolves through time, Lz tends to 
drift, making the meaning of the parameter Q less precise. 
Therefore, we will look for steady solutions with no net ver- 
tical angular momentum relative to the rotating frame. We 
achieve this by calculating the drift in the value of Lz af- 
ter each iteration and then introducing a correction in the 
form of an equivalent rigid body rotation in the opposite 
direction. This alters the evolution of the PDEs (slightly), 
but steady states are correct solutions of the PDEs. In our 
oscillatory solutions we remove the constraint from Lz and 
follow its time evolution. This is discussed in Section [3.91 

A uniform, vertical magnetic field is used as initial con- 
dition. For a nonrotating cylinder (Q = 0) the azimuthal 
magnetic field is perturbed (Section 13. 1|) . For a finite angu- 
lar velocity (11/0) the evolution of the plasma is triggered 
by starting the quiet, nonrotating plasma with a finite Q 
and no plasma perturbation. Both these initialisations en- 
sure that = at the start of the numerical simulations. 

The density does not in principle satisfy a boundary 
condition, but we impose the value of the normal derivative 
of p obtained from the Navier-Stokes equation (g]). 
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Figure 1. The diagnostics used to describe the numerical results: 
1. Potential magnetic field lines; 2. Magnetic field lines; 3. Veloc- 
ity field in the (r, z) plane; 4. Temperature fluctuation relative 
to the unperturbed state; 5. Density contour lines; 6. Magnetic 
fleld direction and strength; 7. Azimuthal current density; 8. Az- 
imuthal velocity fleld; 9. Azimuthal magnetic fleld; 10. Current 
density in the (r, z) plane. All colour scales have red as maximum 
and blue as minimum. The colour green represents zero. 




Figure 2. No rotation with the azimuthal velocity and magnetic 
fleld components approaching zero asymptotically. Consequently 
the right hand side of Figure [T] is not included. The absolute 
temperature variation has a maximum of max IT] = 3.1, while 
the azimuthal current density lies in the interval (-200,350). 



The numerical co de was developed spec i fically for this 
type of calculation (jHurlburt Ruckhdgel l2000h . Sixth- 
order compact finite differencing is used, which reduces at 
the boundaries to fifth-order accuracy for first-order deriva- 
tives and fourth-order accuracy for second-order derivatives. 
The grid intervals were chosen to be equal in the r and z 
directions, with 240 grid points in the horizontal and 80 in 
the vertical for the majority of calculations. The time evo- 
lution obtained fourth-order accuracy through a modified 
(explicit) Bulirsch-Stoer integration scheme, with the time 
step limited by the Courant condition (using the maximum 
sound and Alfven speeds, as well as the thermal diffusive 
limit), multiplied by a safety factor of 0.5. 



3 NUMERICAL RESULTS 

Unless otherwise stated, the results shown here have been 
obtained with the following parameter values: R — 10^, Q — 
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Figure 3. Constant angular velocity Q = 0.1. The time independent solution shows two convection cells in the radial direction instead of 
the one in the result obtained with no rotation (Figure [2]) . The diagnostics are described in Figure [T] with the azimuthal current density 
G (—189, 365) and the current density in the (r, z) plane jr G (—30, 109) and jz G (—120, 106). The temperature variation is such that 
max \T\ = 2.88, and the measured max |iir| = 2.36, max = 0.87, max \uz\ = 2.35, and max l^^l = 3.96. 





radial distanae 



radial distanae 



Figure 4. Radial profile of Ur with fl = 0.1. The solid line is at 
depth 0.25, the dotted line at 0.5, and the dashed line at 0.75. 



Figure 6. Radial profile of the axial or vertical velocity Uz with 
Q = 0.1. The three lines have the same meaning as in Figure 0] 




Figure 5. Radial profile of azimuthal velocity with ^2 = 0.1. 
The lines have the same meaning as in Figure 0] 



Figure 7. Radial profile of azimuthal magnetic field 5^ with 
fl = 0.1. The three lines have the same meaning as in Figure 0] 
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32, cr = 1, Co = 0.2, 6* = 10, m = 1, 7 = 5/3 and T = 3. The 
results are presented in the format given in Figure [1] 



3.1 No rotation 

When no is rotation present = or Ro cxo), the az- 
imuthal magnetic field is perturbed initially and the plasma 
is allowed to evolve through time. The solution reaches the 
state depicted in Figure [21 with the explanation of the diag- 
nostic given by the left hand box in Figure [1] The solution 
described here is typical of the plas ma state when the ax- 
isymmetric cylinder is not rotated ([Hurlburt Rucklidgd 
|200QI : iBotha. Ruckhdge Hurlburtll2006l ) and it provides a 
convenient base against which to compare results obtained 
with rotation. 

From an initial vertical magnetic field, the convection 
sweeps the magnetic field towards the central axis where 
it forms a flux tube. A large anticlockwise convection cell 
forms next to the tube, flowing towards the flux tube at the 
top of the numerical domain. This flow direction keeps the 
magnetic field confined to the central axis. The temperature 
is time dependent in that a cold plasma blob forms at the 
top next to the magnetic flux tube, is convected down the 
side of the tube, only to dissipate as it is convected along the 
bottom boundary. Figure [2] shows a new cold plasma blob 
forming at the top while the remnants of the previous cold 
blob is still visible at the bottom, moving towards the outer 
boundary. This temperature oscillation has a period of ap- 
proximately 1.275 time units, which corresponds roughly to 
half the circulation time around the convection cell. The up- 
per layers of the plasma are heated by the upflow next to the 
outer boundary and the resulting hot plasma blob is time in- 
dependent. The azimuthal current density has its maximum 
value (in both directions) next to the flux tube where the 
magnetic field gradient is the highest. Azimuthal flow and 
magnetic structure are not generated spontaneously. The az- 
imuthal components, generated by the initial perturbation, 
are small and decay exponentially as the solution evolves 
through time. 

Periodic oscillations, an example of which is the time de- 
pendent temperature next to the magnetic flux bundle, are 
famil i ar from Rayleigh-Benard convection (^ Clever Bussd 
Il995l ).[j ones Galloway (1993.) found periodic oscillations 
for a Boussinesq fluid in an axisymmetric cylinder. As must 
be expected, as in our case they found no spontaneous gen- 
eration of azimuthal velocity or magnetic field components. 

Given the model's temperature and density profiles in 
([!)) and ([2]), the sound speed increases and the Alfven speed 
decreases with depth. The inflowing layer at the top of the 
box is deeper than the outflowing layer at the bottom. This 
is ascribed to the fact that the total radial momentum in the 
system is zero. The higher density at the bottom leads to a 
shallower outflowing layer with lower radial velocities trans- 
porting the same momentum outwards as what the deeper 
top layer with higher velocities transports inwards. The sim- 
ulation runs with a maximum Mach number of approxi- 
mately 1. The time step is limited by the thermal diffusivity, 
with the dimensionless thermal conductivity K calculated 
using ()12|) . This constraint on the time step is true for all 
the numerical simulations in this paper. 



3.2 Introduce rotation 

By introducing a constant angular velocity of Q = 0.1 
(Ro = 77.5) we obtain the solution in Figure [3] The one 
convection cell in the case of no rotation (Figure [2|) has split 
into two cells, i.e. a finite Q reduces the characteristic wave- 
length of the convection in the radial direction. This effect o f 
rotation on convection is well known (IChandrasekharlll96lh . 
The convection cell next to the flux tube always has an an- 
ticlockwise flow direction with an inward flow at the top, so 
that it forms a collar which forces the magnetic field toge ther 
at the central axis flBotha. Rucklidge k HurlburtI 120061 ). 

By introducing a finite Q, the centrifugal term in equa- 
tion Q provides a force in the radial direction. This man- 
ifests as a change in density contours, which go from being 
approximately horizontal without rotation (Figure [2]) to be- 
ing slanted at the outer wall (Figure [3]). This boundary effect 
is localized and does not affect the solution in the domain 
interior. The treatment of the outer boundary is discussed 
in more detail in Section [3.81 

There is no time dependence in the solution. The radial 
profiles of the velocities are given in Figures |4] to [6] at three 
depths. All three velocity components are of the same order 
of magnitude. The radial velocity (Figure |4]) shows the two 
cells circulating in opposite directions, as well as the fact 
that the speed is higher in the upper part of the numerical 
domain. 

The azimuthal velocity (Figure [5]) shows that the 
plasma inside the strong magnetic field of the flux tube ro- 
tates as a solid body, with maximum rotation on the outside 
edge of the tube. In the convection area of the solution the 
rotation is in the form of a vortex with the azimuthal velocity 
gradually falling away with radius. This rotation pattern is 
uniform throughout the depth of the box and compares well 
with observations that show t he largest azimutha l velocities 
are located in the penumbra (jBrown et al.ll2003l ). One can 
fit the profile with a Rankine vortex, described by 



— — tor r ^ it, 

-ft 



VoR 



(19) 



for r > R, 



with R the magnetic flux tube radius and Vb = max(i^0). 
An observer in the rotating reference frame of the cylinder 
will measure an azimuthal velocity profile of 



(20) 



The values of Vb and R measured for r2 = 0.1 are presented 
in Table [T] and the radial profile of v^^ in Figure [8] Rank- 
ine vortices are used regularly to model tropical cyclones on 
Earth. Helioseismic measurements of fiow around sunspots 
in the upper convection zone show a strong resemblance to 
the fiow of hurricanes on Earth (Zhao k Kosovichev 20o3). 
It is a happy coincidence t hat Herschel thought of sunspots 
as large cyclonic storms ([Thomas Weisslll992l l. In our 
model the Rankine vortex makes physical sense. Convection 
is suppressed where the magnetic field is strong. The radial 
dependence of the azimuthal velocity in these regions is that 
of a rotating rigid body. Since the region experiencing rigid 
rotation corresponds to the magnetized region in all cases, 
we deduce that magnetic effects are responsible for the rigid 
body rotation. A vortex exists around the fiux tube. Angular 




radial distanae 



Figure 8. Radial profile of a Rankine vortex inside a reference 
frame rotating with ^2 = 0.1. v'^ is described by (|20|) and should 
be compared with Figure [5] The lines have the same meaning as 
in Figure m and the Rankine constants are listed in Table [T] 

momentum mixes in axisymmetric convection, which results 
in the free vortex in the field-free convection cells where con- 
vection is strong. The counter flow near the outer wall is a 
consequence of the treatment of Lz in our solution. Since 
our solution has zero vertical angular momentum relative to 
the rotating reference frame, a significant counter flow has 
to occur at the edge in order to balance the peak flow next 
to the flux tube. 

The vertical velocity (Figure [6]) shows the strong down- 
flow at the outside of the magnetic flux tube and at the outer 
edge of the numerical domain, as well as the upflow between 
the two convection cells. Comparing the two downflows, we 
observe that the downflow at the outer edge is stronger than 
that next to the magnetic flux tube. It is essential to use a 
large enough aspect ratio (F) so that the outer boundary is 
removed from the physics around the magnetic flux tube. A 
F = 3 appears to be a reasonable compromise between this 
and the computational limitations. 

The radial and axial magnetic fleld components are con- 
centrated in the magnetic flux tube, with Bz three times 
larger than Br. Figure [7] shows the radial proflle of B^, the 
size of which is an order of magnitude smaller than Br. B(f) 
is conflned mainly to the inner convection cell next to the 
magnetic flux tube. The current, obtained from the mag- 
netic fleld through equation (|9|, reflects the distribution of 
the magnetic fleld. Its azimuthal component is concentrated 
on the outside of the magnetic flux tube where the radial 
gradient in the magnetic fleld is the largest. The radial and 
vertical components are distributed in and around the inner 
convection cell around the azimuthal magnetic fleld maxima 
(Figure[3|). At the top and bottom boundaries the radial cur- 
rent density has local maxima due to the fact that no current 
flows out of the box. 

3.3 Increase rotation 

The convective Rossby number (Ro) associated with the ro- 
tation around the central axis decreases as Q increases. The 
Ro associated with the circulation around the convection 
cells for the case with Q = 0.1, i.e. Ro = 77.5, compares 
well with that of supergranulation in the Sun. Ro of larger 
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Table 1. The measured constants in the Rankine vortex (I19|) . 



n 


z (measured downward) 


Vo 


R 


0.1 


0.25 


0.75 


0.55 




0.5 


0.69 


0.55 




0.75 


0.57 


0.63 


0.2 


0.25 


1.16 


0.68 




0.5 


1.09 


0.68 




0.75 


0.95 


0.73 


0.3 


0.25 


1.27 


0.85 




0.5 


1.24 


0.84 




0.75 


1.13 


0.88 



Q values correspond to even larger-scale flows. In all cases 
the rotation rate is low enough that it should not signifl- 
cantly change the value of t he critical Rayleigh number for 
the on set of convection (see iBrummelL Hurlburt Toomrd 
(|l996l ll and thus the cases exhibit comparable amplitudes. 

As the magnitude of Q increases, the width of the mag- 
netic flux tube increases. Figure [9] shows a time independent 
solution with Q — 0.3. The magnetic fleld strength inside 
the flux tube decreases with increasing width, allowing weak 
convection to form inside the flux tube itself. Figure [9] shows 
that the upflow in the flux tube heats the plasma in the top 
layers of the tube, while very weak outflow forms along the 
top boundary. Eventually, for Q ^ 0.3, the convection in- 
side the flux tube becomes strong enough to break it into 
concentric rings. 

As rotation increases and with it the width of the mag- 
netic flux tube, the magnetic fleld lines forming the flux tube 
straighten. This causes the azimuthal current density to 
decrease in both positive and negative azimuthal directions. 
The position of stays the same: it flows around the flux 
tube, created by large magnetic fleld gradients there. 

Increasing Q also increases the size of the centrifugal 
force in equation (|4]). For Q = 0.1 the density contours in- 
side the flux tube are approximately horizontal (Figure [3]). 
For Q ^ 0.2 density contours become slanted, due to weak 
convection inside the flux tube as well as the centrifugal 
force acting on the plasma. As a result, there is a slight de- 
pletion of plasma at the top of the magnetic flux tube near 
the central axis, while the area of density variation along the 
inside edge of the flux tube increases. Figure [9] shows that 
inside the convection cells the increase in Q causes a slight 
depression of density contours. 

The plasma rotates as a Rankine vortex for all values 
of Q. Where the magnetic fleld is strong enough to sup- 
press convection, the flow is a forced vortex in the form of 
rigid body rotation, while in the fleld-free convection region 
we observe a free vortex. The 1/r dependence in the convec- 
tion region corresponds to homogeneous angular momentum 
that is caused by the effective mixing by the convection. As 
the width of the magnetic flux tube increases with the in- 
crease in Q, the radius of the forced vortex also increases 
(Table [1]) . The radial proflles of V(f) for different Q values 
can be compared in Figures [5] and 1101 They show that as 
Q increases, the maxima next to the magnetic flux tubes 
increase as well. (See also Table [1]) To maintain the initial 
Lz = 0, the count erflow at t he outer wall increases in sym- 
pathy. This is in contrast to I Jones k Gallowavl (|l993l ). who 
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Figure 9. Constant angular velocity Q = 0.3. The convection inside the magnetic flux tube is stronger and the flux tube wider than 
in th case with Q = 0.1 (Figure [3]). The growth in width is at the expense of the inner convection cell. The weak convection creates a 
temperature signature inside the flux tube. The measured max Iw^l = 1.5, max l^^l = 8.0, max \T\ = 2.9, G (—143, 220), jr G (—61, 75), 
and jz e (-201,111). 



0.5 



0.0 



-0.5 




15 - 



radial distanae 



Figure 10. Radial proflle of azimuthal velocity with Q = 0.3, 
obtained from Figure [9] The lines have the same meaning as in 
Figure |4] 



found a retrograde flow at the central axis for large Q val- 
ues in a Boussinesq fluid. To generate a retrograde flow near 
the central axis, we had to change the temperature bound- 
ary condition at the lower boundary in our model (Section 
13. 6p . The rigid body rotation inside the magnetic flux tube 
is perturbed when the weak convection inside the tube be- 
comes strong enough to influence the local magnetic field. 
Figure [9] shows the strength of the convection inside the flux 
tube, while Figure [10] shows the deviation from rigid body 
rotation. This deviation increases deeper in the numerical 
domain where convection is stronger. 

The azimuthal magnetic field tends to be located in the 
convection cells closest to the central axis. In the case of 
Q — 0.1 this is in the collar flow around the flux bundle 
(Figure [3j). For Q ^ 0.2 a small convection cell starts to 
form inside the flux bundle at its base, due to weak convec- 
tion inside the flux bundle (Figure [9}. As this cell grows in 
strength, the amplitude of located in it grows in strength 



^ 
0.0 



0.1 0. 
Q 



0.3 



Figure 11. The behaviour with increasing rotation is presented 
using the following notation: Br is green plus signs connected by 
a solid line; B^j^ is blue stars connected by a dot-dashed line; Uz is 
cyan diamonds connected by a dotted line; Ur is magenta crosses 
connected by a dashed line; is black triangles connected by a 
triple-dot-dashed line. The peak values are plotted in each case. 



relative to the 5^ in the collar flow outside the flux bundle. 
The directions of in the small cell inside the flux bundle 
and that of in the collar flow are anti-parallel. This corre- 
sponds to the direction of flow of the convection cell. Figure 
[9] shows that a clockwise convection cell has a point- 
ing in the positive direction (i.e. into the page), while i?^ 
in an anticlockwise cell points in the negative (f) direction 
(i.e. out of the page). This is caused by the interaction of 
the magnetic fleld with the velocity in the flrst term on the 
right hand side of equation (|7|). The current surrounding 
the local maxima of i?^ has the same direction as the local 
convection, since it is calculated using equation Q. 

The behaviour when rotation is increased is summa- 
rized in Figure 1111 As Q increases the magnetic flux tube 
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Figure 12. Results with Q = 128 and Q = 0.3. The stronger magnetic field widens the magnetic flux tube at the expense of the size 
of the convection area, as seen when compared to Figure [9] The strength of convection and the size of the azimuthal quantities are 
reduced. Maxl^^l is located in the anticlockwise convection cell next to the flux tube. The measured maxlw^l = 1.0, maxl^^l = 2.9, 
max |T| = 2.7 and G (—65, 125). The range of the current density in the (r, z) plane is jV G (—13, 54) and jz G (—74, 33). 




radial distanae 



Figure 13. Radial proflle of corresponding to Figure [T2] with 
Q = 0.3 and Q = 128. The lines have the same meaning as in 
Figure |4] 



widens and the field lines straighten and become more ver- 
tical, which leads to a decrease in the radial component of 
the magnetic field. At the same time the azimuthal veloc- 
ity and magnetic field components increase, being driven by 
Q. Compared to these changes, the radial and axial velocity 
components stay relatively stable. All velocity components 
increase in absolute value as Q increases. 



3.4 Increase magnetic field strength 

From previous numerical studies it is known that an increase 
in the magnetic field increases the width of the magnetic flux 
tube for nonrotating solutions (Hurlburt & Rucklidge 2000). 
This is also true when rotation is present, which can be seen 
when Figure [9] with Q = 0.3 and Q = 32 is compared with 
Figure [12] for Q = 0.3 and Q = 128. The solution is time 
independent and the growth in flux tube width takes place 



at the expense of the radial dimensions of the convection 
cells. 

The magnetic flux tube retains its configuration, with 
an anticlockwise convection cell holding the flux tube in 
place. The stronger magnetic field suppresses the weak con- 
vection inside the magnetic flux tube, which was present 
when Q = 32 (Figure [9]). 

As the area with strong magnetic field becomes wider, 
the field- free convective region is compressed. The decrease 
in the field-free area is accompanied by lower flow veloci- 
ties of convection. Here the maximum Mach number of the 
solution is 0.8, while max (Mach) = 0.9 for Q = 32. The 
maximum measured azimuthal velocity for Q = 128 (Figure 
fT2l) is also 2/3 of what it is for Q = 32 (Figure [9]). 

The weaker flow in the convection cell around the mag- 
netic flux tube means the field lines are less compressed when 
compared to the case when Q = 32 (Figure [9]). This leads 
to lower gradients at the flux tube's edge, which in turn im- 
plies a lower azimuthal current density flowing around the 
flux tube, since jcf, is calculated using equation (|9]). 

The azimuthal flow of this solution fits that of a Rankine 
vortex. By suppressing the weak convection inside the flux 
tube that is present for Q = 32, the plasma flow inside the 
flux tube becomes more like rigid body rotation. (Compare 
Figures [10] and [13]) The maximum U(j) next to the flux bundle 
is lower for Q = 128 due to the lower levels of convection 
in the solution. This is also true for the counter flow at the 
outer wall. 

The azimuthal magnetic field has its maximum in the 
convection cells closest to the central axis. In Figure [9] with 
Q = 32, there exists a small clockwise cell at the base of the 
flux tube. This cell is strong enough to contain a significant 
part of B(f,, with the anticlockwise collar flow containing an 
anti-parallel component. Here, for Q = 128 (Figure [12]), 
the small clockwise cell inside the flux tube is suppressed, so 
that max \B(f)\ is mainly located in the anticlockwise collar 
flow. 
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radial distanae 

Figure 15. Radial profile of li^ corresponding to Figure [M] with 
fl = 0.1, cr = 0.3 and Q = 32. The lines have the same meaning 
as in Figure |4] 



3.5 Lower Prandtl number 

Decreasing the Prandtl number a causes the convection de- 
scribed by the steady solution to become more vigorous, as 
can be seen when Figure [3] (with a = 1) is compared to 
Figure [Ml (with a — 0.3). For a — 1 the maximum Mach 
number in the solution is 0.9, while for a = 0.3 we have 
max (Mach) = 1.7. The lower a value brings the simulation 
closer to the physical conditions in the upper convection 
zone. The dimensionless thermal conductivity defined 
by (O, changes from 4.9 x 10~^ for a = 1 to 1.5 x 10~^ 
for cr = 0.1. For a = 0.3 we have i^^ = 8.9 x 10~^ and with 
= 0.1 we have a convective Rossby number of Ro = 42.4. 
For lower Prandtl numbers the inner convection cell in- 
creases its size at the expense of the width of the magnetic 
flux bundle, and to a lesser degree the width of the convec- 
tion cell next to the outer boundary. Not only is the inner 
cell larger, but the velocity amplitudes have higher maxi- 
mum values. However, the relative differences between the 



velocity components are independent of the value of cr, with 
the azimuthal component approximately a third the size of 
the radial and axial components. The azimuthal flow takes 
the form of a Rankine vortex (Figure [15]) • The rigid body 
rotation inside the flux bundle is faster than when cr = 1 
(Figure [5|), with a higher maximum next to the bundle. The 
flux bundle is also narrower, which means that to maintain 
the initial = during the simulation, the counter flow 
at the outer boundary needs to be only slightly larger than 
when cr = 1. 

The stronger convection pushes the magnetic flux into a 
thinner flux bundle at the central axis, so that the strength 
of Bz increases for lower values of cr. The size of relative 
to Br stays approximately the same for all values of cr. As in 
the case for a — 1 (Figure[3|), the azimuthal magnetic fleld is 
mostly located in the inner convection cell, with its maxima 
next to the magnetic flux bundle. The stronger convection 
also causes the curvature and gradients of the magnetic field 
lines in the (r, z) plane to increase, which increases the size 
of the azimuthal current density obtained from equation (|9|) . 
The position of max Ij^I stays the same. 

The effects of the enhanced convection are visible in 
the density contours. For cr = 1 the density contours in 
the convection cells are approximately horizontal (Figure [3|) 
while for cr = 0.3 significant variation in the radial direction 
occurs f Figure [T^ . At the central axis, where the convection 
is suppressed by the strong magnetic field, the contour lines 
stay approximately as they were with cr = 1 and its lower 
convection strengths. 

There are changes in the temperature profile of the so- 
lution. The stronger upflow between the two convection cells 
leads to a larger variation from the original heat proflle in 
the upper layers of the solution, as can be seen when max |T| 
in Figures [3l and [T^ are compared. Also, for lower cr values 
the thermal diffusion rate becomes more signiflcant, which 
reduces the radial extent of the heated plasma above the 
upflow. 
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Figure 16. Results with a constant dT/dz used as bottom boundary condition. The parameters are ^2 = 0.1, Q = 32 and cr = 0.1. Weak 
convection forms inside the wide magnetic flux tube. The overall convection is lower and the inner convection cell smaller compared to 
solutions using a constant T lower boundary condition. The measured max \U(j)\ = 0.8, max = 3.2, max \T\ = 1.4 and j^^ G (—73, 105). 
The current density in the (r, z) plane is jV G (—27, 23) and jz G (—53, 50). 



Table 2. Changing the lower temperature boundary condition with Q = 32, m = 1, F = 3, ^ = 10. 



n 


cr 


Lower boundary: T constant 


Lower boundary: dT/dz constant 






max(Mach) 


T dit z = l 


max(Mach) 


T range at z = 1 





1.0 


1.2 


11 


0.6 


(9.34,10.19) 


0.1 


1.0 


0.9 


11 


0.6 


(9.33,10.42) 


0.1 


0.3 


1.7 


11 


1.1 


(9.05,10.43) 


0.1 


0.1 






1.3 


(8.59,10.20) 



3.6 Temperature prescription at bottom 
boundary 

To determine the extent to which the bottom boundary is 
influencing the result, we changed the temperature prescrip- 
tion on this boundary from a constant value T to a constant 
dT/dz. From equation ([1]) we set dT/dz = ^, so that the 
heat flux K6 stays equivalent to the heat flux for a con- 
stant T boundary co ndition. A linear stability analysis by 
iHurlburt, Toomre Massaguerl (|l984[ ) determined that this 
change in boundary condition results in halving the criti- 
cal Rayleigh number for the onset of convection, and hence 
one would expect somewhat more vigorous convection for 
the same Rayleigh number, all other aspects of the solu- 
tion being equal. However, the numerical results discussed 
in this section show that for this highly nonlinear system, 
the change in the lower boundary condition leads to lower 
convection levels. 

For low Prandtl numbers, such as a — 0.1 in Figure 
[T6l the basic configuration of two convection cells and a 
central magnetic flux tube remains as before. The radial 
profile of the azimuthal flow is that of a Rankine vortex. 
The strength of convection outside the magnetic flux tube 
is lower than for the bottom boundary condition of constant 
temperature, while the magnetic flux tube is wider with very 
weak convection inside it. 

For Prandtl numbers of a ^ 0.3 the solution changes 
into one convection cell that is outflowing at the top. This 



new flow direction does not provide an eflRcient collar to con - 
tain the magnetic flux teotha, Rucklidge Hurlburtll2006l ) , 
so that the magnetic field spreads out in the radial direction 
rather than being contained at the central axis. Under these 
circumstances it is possible to have horizontal magnetic field 
lines above the main convection cell. Figure [17] gives the so- 
lution for (7 = 1, while Figures [18] to [20] give the radial 
profiles of for various parameter values. 

Figure [T71 shows the solution with cr = 1. A solution in 
the same parameter space but with constant temperature at 
the bottom boundary is given by Figure [3] When compar- 
ing the two solutions, it is clear that the level of convection 
is lower in the case of constant dT/dz boundary condition. 
This is shown explicitly in Table [2] the maximum Mach num- 
ber is lower when a constant dT/dz is used. It also shows 
in the fact that the maximum azimuthal velocity is weaker 
in Figure [17] than in Figure [3] This is true for all choices of 
parameter values. The maximum Mach number in the solu- 
tion increases as a decreases, in line with the discussion in 
Section 13.51 and the fact that more heat fiows through the 
system. In Table [2] the solution for a = 0.1 and constant 
dT/dz does not have a counterpart when a constant tem- 
perature boundary condition is used, because the convection 
becomes too vigorous and large shocks form that terminate 
the numerical simulation. This agrees with the conclusion 
that a constant dT/dz at the lower boundary leads to lower 
convection levels. 

The weaker convection allows the magnetic fiux tube to 
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Figure 17. Results with constant dT/dz bottom boundary condition and the same parameter values as Figure [T6] but with (7 = 1. One 
convection cell forms with outflow along the top boundary. This allows the magnetic flux to spread radially, with weak convection forming 
near the central axis inside the flux tube. The measured max|wr| = 0.8, max|W(^| = 0.5, maxlw^^l = 0.8, maxl^^l = 3.7, max|T| = 1.1 
and G (-59, 149). The current density in the (r, z) plane is jV G (-20, 22) and jz G (-134, 50). 




Figure 18. Radial proflle of with ^2 = 0.1, a = 0.3, Q = 32 
and a constant dT/dz as bottom boundary. The solution has a 
conflguration similar to Figure [TTI but the weak convection inside 
the flux tube has higher amplitudes. These are responsible for the 
large perturbations in the interval ^ r ^ 1. The lines have the 
same meaning as in Figure 0] 



be wider, so that the value of max \ Bz\ is lower. The differ- 
ent flow pattern when a ^ 0.3 also contribute to the lower 
vertical magnetic component, in that the magnetic field now 
has a significant horizontal component (Figure [TT]). Another 
consequence of the weaker convection is lower magnetic field 
gradients in the (r, z) plane. This leads to lower levels of az- 
imuthal current density, which can be seen when Figures [3] 
and [TTI are compared. In the case for a = 0.1 (Figure [T6|) the 
peak current density is positioned around the magnetic flux 
tube close to the midplane, where it is for all results with a 
collar flow around the flux tube, while for a ^ 0.3 (Figure 
I17p the maximum current density is located at the base of 
the flux tube, where magnetic field gradients are largest. 

As in the case for a constant T lower boundary, the 
plasma rotates as a rigid body (or forced vortex) where the 



Figure 19. Radial proflle of corresponding to Figure fTTI with 
Q = 0.1, (J = 1, and Q = 32. Figure |4] deflnes the line notation. 
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Figure 20. Radial proflle of with Q = 0.1, a = 0.3, Q = 256 
and constant dT/dz. The solution has a conflguration similar to 
Figure [TT] The lines have the same meaning as in Figure 0] 
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magnetic field is strongest, while a free vortex forms in the 
convection area where the magnetic field is weaker. For a ^ 
0.3 the central magnetic flux tube is still present, with an 
additional component that stretches horizontally above the 
dominant convection cell. This means that at the bottom of 
the numerical box, where one has a well defined flux tube 
(Figure [T7|) . the radial profile of the azimuthal flow looks 
most like a Rankine vortex. (See Figures [18] for a = 0.3 and 
[19] for a = 1.) Moving higher up in the numerical domain, 
the width of the magnetic flux tube increases and with it 
the radius of the forced vortex, with the free vortex in the 
convection area occupying less space. At the top of the box 
the magnetic field influences the azimuthal flow so much 
that most of the free vortex flow is distorted. The strong 
convection associated with low a values fSection l3.5|) enables 
weak convection to form inside the magnetic flux tube. This 
weak flow serves as a perturbation to the rigid body rotation 
inside the magnetic flux tube around the axis (Figure [18] for 
a = 0.3). 

By increasing a one weakens the convection in the so- 
lution, as discussed in Section [3.51 and reduces the heat flux 
through the system. This allows the magnetic field to be- 
come more uniform along the central axis with less pertur- 
bations in the rigid body rotation that occurs there. Figure 
1171 shows an example of this for a = 1. An interesting phe- 
nomenon occurs when cr = 1. In this case the plasma inside 
the magnetic flux tube slows down to almost zero (Figure 
I19p . The plasma in the field- free convection area still forms 
a free vortex, which means the azimuthal flow grows from 
a small value to its maximum next to the flux tube over 
a small distance, before it tails off into the usual free vor- 
tex profile. The maximum flow in the vortex occurs close 
to the base of the flux tube, where the convection area has 
the lowest level of magnetic field. The fact that the slowly 
flowing plasma inside the flux tube is still in the same di- 
rection as the free vortex is a coincidence. The plasma flow 
inside the magnetic field area is highly sensitive to the val- 
ues of Q and a. Some values give a retrograde rotation at 
the central axis, in the same direction as the counter flow 
near the outer wall. Other values give a prograde rotation in 
the direction of the free vortex, but with huge perturbations 
in the rigid body rotation, while others will give a flux tube 
that rotates partly prograde and partly retrograde. As one 
example of what can occur, we plot the radial profile of 
for the values Q = 256 and a — 0.3 in Figure 1201 Here the 
plasma in the whole of the flux tube rotates retrograde, as 
well as most of the plasma in the magnetic flux layer on top 
of the one large convection cell. (The solution has the same 
configuration as Figure [TT] only with the width of the flux 
tube wider due to the higher Q value.) Notice in Figure [20] 
that the free vortex still rotates prograde with a counter flow 
next to the outer wall, similar to Figure [19] when Q = 32 
and cr = 1. 

For a constant T at the lower boundary, the azimuthal 
magnetic field formed in the inner convection cell closest 
to the central axis. Here, for a = 0.1 (Figure [T6|) . weak 
convection inside the magnetic flux tube allows a small cell 
to form at the base of the flux tube. This cell is closest 
to the central axis and carries a large part of B(f) in it. The 
much stronger convection cell forming the collar flow around 
the magnetic flux tube contains the rest of B(j). These two 
cells have opposite meridional circulations, so that the 



components in them are anti-parallel to each other. For the 
cases where a ^ 0.3 (Figure [TT]), only one cell dominate 
the convection area. The azimuthal magnetic field is located 
inside this cell, but with its maximum value next to the 
magnetic flux tube. It is interesting to note that max(B(/,) 
is situated towards the top of the numerical domain, while 
it is towards the bottom of the domain for constant T lower 
boundaries. This corresponds to the direction of flow in the 
convection cell containing in each case. 

The influence of the outer wall is discernible with con- 
stant dT /dz at the lower boundary, when only one clockwise 
convection cell forms in the solution (i.e. for a ^ 0.3). At the 
outer wall a local max(B0) forms, generating its own cur- 
rent around it in the {r^z) plane (Figure [TT]). The heat flux 
through the system with cr < 1 is higher than for a — 1 and 
the convection stronger, so that the magnetic field is less able 
to concentrate next to the outer wall. One can also see the 
influence of the outside wall in the azimuthal velocity pro- 
flle; the amplitude of the counter flow next to the outer wall 
diminishes sharply at the wall. This effect becomes larger 
as the size of the local max(B<^) at the outer wall increases. 
(Compare Figures [T8] and [T9] where a increases, as well as [18] 
and [20] where Q increases.) In all our results this boundary 
effect is highly localized and does not influence the solution 
deeper in the numerical domain. In Section 13.81 the treat- 
ment of the outer boundary is discussed. 

Figures [21] and [22] show that the temperatures at the 
top and bottom domain boundaries are lower when con- 
stant dT /dz is used, compared to a constant T at the lower 
boundary. Inside the magnetic flux tube the temperature 
near the top of the domain (Figure [2T]) is largely indepen- 
dent of the bottom boundary condition. Where the convec- 
tion dominates outside the magnetic flux tube, the temper- 
ature variation is lower than for a constant T at the bot- 
tom. Figure [22] shows the temperature variation at the lower 
boundary. Similar to the top of the domain, the tempera- 
ture inside the magnetic flux tube is least affected by the 
boundary condition, although the temperature is lower along 
the whole radial length for dT/dz constant. Figure [22] and 
Table [2] show that the variation of the temperature along 
the bottom boundary is substantial. This observation may 
ex plain the difference between the linear stability results 
by iHurlburt, Toomre &: Massaguerl (p.984) and the nonlin- 
ear behaviour, as mentioned in the beginning of this section. 
Although the heat flux into the domain stays the same, a 
lower temperature along the bottom boundary will drive 
the convection less vigorously. To investigate how sensitive 
the velocity amplitudes are to the lower bottom tempera- 
ture, we increased the bottom boundary temperature and 
the heat flux through it by increasing 0. This, however, has 
repercussions throughout the system, as discussed in Section 
[32] 

3.7 Increase stratification in numerical domain 

Increasing the value of from 10 to 20 is felt through- 
out the system. From equations ([1]) and ([2]) we see that 
the stratiflcation doubles. Through equation ()12|) the value 
of the dimensionless thermal conductivity K changes from 
4.9 X 10"^ for 6> = 10 to 1.3 x 10"^ for = 20. This means 
the heat flux through the system [KO) is raised flve fold, 
as mentioned in Section 13.61 It also implies that the con- 
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Table 3. Changing the stratification of the domain (0) with cr = 1, m = 1, F = 3. 



n 


Q 


Lower boundary 


= 


10 


= 


20 








max(Mach) 


T at z = l 


max(Mach) 


T at z = 1 


0.1 


32 


dT/dz constant 


0.6 


(9.33,10.42) 


0.8 


(17.81,19.53) 


0.3 


128 


T constant 


0.8 


11 


1.0 


21 



1 .50 r 

1 .40 h 
1 .30 r 
1 .20 h 
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Figure 21. Top boundary temperature with constant dT/dz at 
bottom boundary. The soHd hue is cr = 0.1 (Figure[T6]), the broken 
hue (7 = 0.3 and the dot-dashed hue a = 1 fFigure fTT)) . The dotted 
hne is a = 0.3 with a constant bottom temperature (Figure 
added as reference. 




radial distance 



Figure 22. Temperature at bottom boundary with dT/dz con- 
stant. The hues have the same meaning as in Figure [2T1 The extra 
(triple-dot-dashed) hne corresponds to cr = 1 and ^ = 20 (Figure 
I24|) . as discussed in Section [3771 To fit into the graph, 10 has been 
subtracted from this temperature. 



vective Rossby number ()13p increases one order of magni- 
tude while Q stays unchanged. Table [3] shows that the max- 
imum Mach number in the solution increases for both tem- 
perature prescriptions. The linear stability properties are 
changed both by increasing the critical Rayleigh number for 
the onset of convection, due to the ch ange in stratification 
(|Hurlburt, Toomre Massaguerl 1 1984 ), and by incr easing 



the mean value of in the domain (jWeiss et al]ll99d l. 

Increasing does not change the basic configuration 



of the numerical solution. The magnetic flux tube forming 
at the central axis remains intact, as well as the convec- 
tion cells in the field-free region. This is true for a constant 
temperature as well as a constant dT/dz bottom boundary 
condition. Increasing increases the strength of convection 
in the (r, z) plane, as measured by the Mach number in Ta- 
ble [3l as well as the width of the magnetic flux tube. The 
latter can be observed in Figure [23] for a constant T at the 
bottom boundary and in Figure [24] for dT/dz constant. The 
wider flux tube leads to lower gradients in the magnetic field, 
which means the azimuthal current density around the tube, 
calculated using equation ([9]), becomes weaker. All the other 
azimuthal quantities (B^ and Ucf,) are also weaker when com- 
pared to results with ^ = 10. 

Figure [23] shows a solution with Q = 128 and constant 
T at the lower boundary. When this is compared to a so- 
lution with the same parameter values and boundary con- 
ditions, but with ^ = 10 (Figure [12]), one observes that the 
wider magnetic flux tube allows weak convection cells to 
form inside it. This convection is strong enough to perturb 
the temperature inside the flux tube in the top half of the 
numerical domain. A careful inspection of the top bound- 
ary shows that these convection cells cause flow along the 
boundary, so that concentric rings start to appear at the 
top. This is a consequence of the axisymmetry in our model, 
as one would expect cellular convection to form inside the 
flux tube. These flows of concentric rings around the central 
axis grow in size as the magnetic flux tube becomes wider, 
which occurs for higher values of Q. The azimuthal magnetic 
field has its maximum value in the strong collar flow next to 
the flux bundle. However, the weak convection cells inside 
the flux tube are defined well enough for B(f) to have sig- 
nificant components inside them (Figure [23]) . The direction 
of each convection cell determines the direction of the local 
B(f, inside it: clockwise convection contains a local maximum 
and anticlockwise a local minimum. The azimuthal velocity 
forms a Rankine vortex, as shown in Figure 1251 The weak 
convection inside the magnetic flux tube perturbs the rigid 
body rotation of the plasma inside the tube. The maximum 
value of U(j) is found next to the outer edge of the flux tube 
and a free vortex forms in the convection area. At the outer 
wall a count erflow forms, as in the case when ^ = 10. A com- 
parison between the two sets of results ( Figure [13] for ^ = 10 
and Figure [25] for — 20) shows that max(i^0) is lower and 
situated farther from the central axis for ^ = 20. This means 
the counter flow at the outer boundary, generated because 
our solution has zero vertical angular momentum relative 
to the rotating reference frame, is approximately the same 
strength for both values. 

Figure [23] shows a solution with Q = 32 and dT/dz 
constant at the lower boundary. This should be compared 
to Figure [17] that has the same parameter values and bound- 
ary conditions, but with ^ = 10. This comparison shows that 
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Figure 23. Results with Q = 128, Q = 0.3 and = 20 and a constant T as bottom boundary. Compared to the results for ^ = 10 (Figure 
I12|) , the magnetic field is more radially dispersed, allowing convection cells to form throughout the radial domain. As one moves towards 
the outer boundary, the upflows between cells grow stronger, accompanied by stronger heating above them. The measured max = 0.7, 
max \B^\ = 1.3, max \f\ = 4.3 and G (-39, 57). In the (r, z) plane we have > G (-4, 15) and jz G (-25, 12). 




Figure 24. Results with Q = 32, Q = 0.1 and ^ = 20 and constant dT/dz as bottom boundary. Compared to the results for ^ = 10 
(Figure El 5 the magnetic field is more radially dispersed, allowing weak convection cells to form inside the flux tube. The strong downflow 
at the outer edge cools the plasma in the lower layers of the domain. The measured max|wr| = 1-0, maxl^^l = 0.5, maxlw^l = 1-1, 
max \B^\ = 2.0, max \f\ = 2.2, G (-41, 74), > G (-10, 7) and jz G (-41, 21). 



the magnetic flux tube is wider for ^ = 20 and that weak 
convection occurs inside the tube. This convection is not 
strong enough to significantly heat the upper layers of the 
numerical domain. In fact, the strong downflows in the so- 
lution cools the lower layers of the numerical domain much 
more than when ^ = 10 (Figure [22]) . The azimuthal mag- 
netic field is located inside the large convection cell next to 
the magnetic flux tube. As in the case for 6 — 10, the maxi- 
mum value of B(f) is located close to the flux tube. Inside the 
field-free convection areas we observe a free vortex, with its 
maximum next to the edge of the flux tube and a counter 
flow next to the outer wall (Figure [26]) . The rotation of the 
plasma inside the flux tube is that of a rigid body, but in 
the opposite direction from the direction of the free vortex 
around the tube. We also see that the plasma inside the hor- 



izontal magnetic field on top of the convection zone rotates 
retrograde, i.e. in the same direction as the counter flow at 
the outside wall. This flow pattern is not a surprise, given 
the fact that we could generate retro flows at the central 
axis and the top of the numerical domain by playing with 
the parameter values in the set of results with ^ = 10. (See 
the discussion in Section [3.61 ) 

The variation in temperature is much larger for these 
cases with ^ = 20 than in the comparable cases with 6 = 10. 
For a constant T at the bottom boundary (Figure [23]) the 
heating occurring at the top of the numerical domain due to 
the strong upflow between the two large convection cells is 
three times larger than the heating for ^ = 10 (Figure [12]). 
In contrast, the result with dT/dz constant at the bottom 
boundary (Figure [2^ shows that the strong downflow next 



16 Botha, Busse, Hurlhurt & Rucklidge 
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Figure 25. Radial profile of corresponding to Figure [23] with 
Q = 0.3, Q = 128 and 6 = 20. This should be compared with the 
case when ^ = 10 in Figure [13] The lines have the same meaning 
as in Figure [4] 
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Figure 26. Radial profile of corresponding to Figure [24] with 
r2 = 0.1, Q = 32 and = 20. This should be compared to the 
case when ^ = 10 in Figure [TO] The lines have the same meaning 
as in Figure [4] 



to the outer wall cools the lower part of the numerical do- 
main. The amount of cooling is double that which occurs for 
= 10 (Figure [HI). 



3.8 Boundary layer at tlie outer boundary 

When comparing results with a rotating cylinder, one ob- 
serves a collar flow next to the flux bundle and a clockwise 
convection cell at the outer boundary. The size of the cell 
at the outer boundary seems to be robust for the various 
parameter values. The only exceptions are for dT/dz con- 
stant at the bottom boundary and with a ^ 0.3 (Figure 
[T7|) , when the whole convection pattern change. Through- 
out the simulations we have taken care that the convection 
cell at the outer boundary does not influence the physics 
near the central axis and the flux tube. Another effect of 
finite Q is that the density contours become slanted at the 
outer boundary due to the centrifugal term in the Navier- 
Stokes equation (|4]). This effect is much more noticeable for 



a constant T bottom boundary condition than when dT/dz 
is used. In contrast, the azimuthal velocity shows a sharp 
decrease in its value at the outside wall when a dT/dz bot- 
tom boundary condition is used (Figures [TH] to [20]) while the 
outside wall hardly register in the U(j) profile for a constant 
T bottom boundary (Figure llOp . The slanted p contours 
and the decrease in amplitude show how the influence of 
the outer boundary on the solution increases as Q increases. 
Due to the formulation of the problem this is unavoidable, 
but it does not pose a problem as long as these effects stay 
localized at the outer boundary. 

In order to minimize the effect of the outer boundary on 
the solution, it was treated throughout as a slippery bound- 
ary, so that the condition on U(f) is given by (|16p , obtained 
from the off-diagonal elements of the rate of strain tensor 
(pTj) . In this paper the boundary conditions at the outer 
wall were chosen so that the coupling between the numerical 
domain and its outside surroundings is kept to a minimum. 
With boundary conditions (p^ only a vertical current exists 
and the Lorentz force is zero at the outer wall. To measure 
the influence of the outer wall on the solution, we changed 
its magnetic boundary condition to that of a perfect con- 
ductor. In this case no currents exist parallel to the wall, 
with a radial current moving through the outer wall. The 
condition that no vertical current exists leads to 



dB^ 
dr 



B(j) 
r 



(21) 



while the radial magnetic field component stays zero, as in 
([16]). For a perfect conductor the Lorentz force has compo- 
nents parallel to the outer wall, but there is no force across 
the wall. This implies that there is a torque at the outer 
boundary, leading to a contribution to the angular momen- 
tum. 

Changing the outer boundary conditions on B(f) to (|2ip 
changes the solution slightly, but only when the azimuthal 
amplitudes at the boundary become significant when com- 
pared to the solution near the central axis, as was the case for 
dT/dz constant at the bottom boundary and a ^ 0.3 (Fig- 
ure [T7|). When the electrically insulating outer wall (Figure 
I17p is compared with a perfectly conducting wall, the solu- 
tion is only slightly perturbed close to the outer boundary. 
The boundary conditions of the bottom boundary, described 
in (|17p , allow currents parallel to the lower boundary but not 
through it. Ditto for the Lorentz force. As a result, in all sim- 
ulations with a perfect conductor at the outer boundary, the 
largest current entering the numerical domain is situated in 
the bottom right hand corner. We attempted to change the 
magnetic boundary condition on the lower boundary, but 
found that the solution is highly sensitive to any changes 
and becomes numerically unstable. Changing the boundary 
condition of on the outer boundary left the azimuthal 
velocity field intact. 

The difference between boundary conditions (p^ and 
(|2ip was thoroughly tested. We started numerical runs from 
a uniform magnetic field for both sets of conditions, which 
lead to almost identical time independent solutions. The nu- 
merical results presented in this paper were started with (|2ip 
and then continued with (|16p . In all cases no significant dif- 
ference between the numerical solutions could be observed. 
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Table 4. Survey of numerical solutions obtained with T constant at the lower boundary and with parameters R = 10^, Co = 0.2, 
171 = 1, 7 = 5/3, r = 3. The star superscript in the Q column indicates time dependent solutions. 



Q 


Q 


a 





Ro 


max |T| 


max(Mach) 


max l^^l 


max 1 Br 1 


max 1 B(f) 1 


max 1 Bz 1 


jrh range 


0* 


32 


1 


10 


oo 


(2.9,3.1) 


(1.2,1.3) 





(13,16.7) 





(30.5,38.8) 


(-465,555) 


0* 


128 


1 


10 


oo 


(2.8,2.9) 


(1.0,1.2) 





(7.3,7.9) 





(14.6,16.9) 


(-184,205) 





256 


1 


10 


oo 


2.7 


0.9 





5.3 





10.8 


(-80,135) 


0.1 


32 


0.3 


10 


42.2 


3.7 


1.7 


1.6 


17.8 


5.5 


43.2 


(-229,445) 


0.1 


32 


0.6 


10 


59.8 


3.2 


1.1 


1.1 


15.2 


4.5 


37.0 


(-208,393) 


0.1 


32 


1 


10 


77.5 


2.9 


0.9 


0.9 


13.3 


4.0 


31.7 


(-189,365) 


0.1 


128 


0.3 


10 


42.2 


3.3 


1.7 


0.8 


8.6 


1.7 


18.9 


(-107,192) 


0.1 


128 


0.6 


10 


59.8 


3.0 


1.1 


0.6 


7.6 


1.6 


15.9 


(-86,171) 


0.1 


128 


1 


10 


77.5 


2.8 


0.8 


0.5 


6.9 


1.5 


14.2 


(-90,178) 


0.1 


128 


1 


20 


205.5 


4.3 


1.0 


0.3 


3.8 


0.6 


7.9 


(-38,67) 


0.1 


256 


0.3 


10 


42.2 


2.6 


1.1 


0.6 


5.2 


0.6 


11.2 


(-65,103) 


0.1 


256 


0.6 


10 


59.8 


2.8 


1.0 


0.4 


5.0 


1.0 


9.5 


(-69,112) 


0.1 


256 


1 


10 


77.5 


2.7 


0.7 


0.3 


4.7 


0.9 


8.9 


(-61,107) 


0.1* 


256 


1 


20 


205.8 


(3.8,3.9) 


0.6 


(0.1,1.5) 


(2.0,2.0) 


0.2 


(7.0,7.1) 


(-33,39) 


0.2 


32 


1 


10 


38.7 


2.9 


0.9 


1.3 


10.7 


4.2 


24.4 


(-148,271) 


0.2 


32 


1 


20 


102.8 


4.7 


1.2 


1.3 


7.4 


2.7 


17.7 


(-74,153) 


0.2 


128 


1 


10 


38.7 


2.8 


0.8 


0.8 


6.3 


2.6 


12.8 


(-73,149) 


0.2 


128 


1 


20 


102.8 


4.3 


1.0 


0.6 


3.7 


1.1 


7.5 


(-37,64) 


0.2 


256 


1 


10 


38.7 


2.7 


0.7 


0.6 


4.4 


1.7 


8.4 


(-56,101) 


0.2* 


256 


1 


20 


102.8 


(3.7,3.9) 


(0.6,0.7) 


0.3 


(1.8,2.0) 


(0.3,0.4) 


(4.5,4.6) 


(-32,37) 


0.3 


32 


1 


10 


25.8 


2.9 


0.9 


1.5 


8.5 


8.0 


20.3 


(-143,220) 


0.3 


32 


1 


20 


68.5 


4.6 


1.2 


1.5 


6.3 


2.8 


14.3 


(-66,123) 


0.3 


128 


1 


10 


25.8 


2.7 


0.8 


1.0 


5.3 


2.9 


10.6 


(-65,125) 


0.3 


128 


1 


20 


68.5 


4.3 


1.0 


0.7 


3.3 


1.3 


7.8 


(-39,57) 


0.3 


256 


1 


10 


25.8 


2.6 


0.6 


0.7 


3.8 


2.1 


7.2 


(-55,89) 


0.3* 


256 


1 


20 


68.5 


3.7 


0.5 


(1.6,1.7) 


(1.6,1.7) 


0.4 


3.8 


(-31,37) 



3.9 Time dependence 

Increasing the angular velocity Q increases the width of the 
magnetic flux tube and allows weak convection cells to form 
inside the tube, similar to those formed in Figure [23] If the 
value of Q becomes large enough, these cells undergo peri- 
odic motion. For a constant T bottom boundary condition, 
Q = 256 and ^ = 20, the weak convection cells inside the 
flux tube oscillate in the radial direction. As Q increases 
from 0.1 to 0.3, the amplitude of this oscillation increases as 
well. In contrast, for a bottom boundary condition of con- 
stant dT/dz with Q = 0.3 and ^ = 10, a hot blob forms due 
to weak upflows inside the tube. This blob then moves to- 
wards the central axis where it dissipates. A new blob forms 
inside the flux tube and the process repeats itself. This hap- 
pens for Q values of 128 and 256 with Q = 0.3, as well as 
for Q = 256 and Q = 0.2. When the value of is doubled to 
20, the solutions become time independent again. 

When the forced conservation of Lz is lifted, its value 
tends to drift, as discussed in Section 12.21 In the case of 
a time dependent solution, Lz oscillates in sympathy with 
the oscillation in the convection, in addition to its drift. The 
amplitude of the oscillation can thus be expressed in terms of 
an equivalent solid body rotation, which is of order O(10~^). 
This compares to a drift in the value of Lz of order O(10~^) 
per unit time. 



4 SUMMARY 

We have investigated magnetoconvection around a magnetic 
flux bundle in a cylinder, when the cylinder is rotated at a 



constant angular velocity Q. The model uses a compressible 
plasma with density and temperature gradients simulating 
the upper solar convection zone. All the numerical solutions 
that we obtained are presented in Tables |4] and [5] Through- 
out the calculations the maximum velocities are in the (r, z) 
plane, so the the maximum Mach number in these tables are 
a good proxy for Ur and Uz. For time dependent solutions 
we present the range in which the different diagnostics lie. 

With no rotation (Q = 0) and a constant temperature 
at the lower boundary, the solution is in the form of a flux 
tube situated at the central axis, surrounded by a field-free 
annular convection ring that forms a collar around the flux 
tube (Section I3.ip . This magnetic configuration lends itself 
to the description of idealized pores and sunspots. The col- 
lar flow has been meas ured in the convection a r ound both 
phenomena. (See >Botha, Rucklidge HurlburtI (|2006[ ) and 
references in it.) 

The introduction of a constant angular velocity Q 
widens the magnetic flux tube (Section I3.2|) . Other ways 
to increase the tube width are to increase the magnetic field 
strength (Section 13. 4p and to increase the heat flux into the 
numerical domain from below (Section 13. 6p . If the magnetic 
field strength (i.e. Q) is kept constant and the tube width is 
increased by means of one of the above, then the amplitude 
of the vertical magnetic field component in the flux tube 
is lowered. This allows weak convection cells to form inside 
the tube. As Q increases, the flux tube widens and the weak 
convection becomes stronger so that eventually concentric 
rings appear at the top of the numerical domain (Figures 
[9] and I23p . In a fully 3D model one would expect cellular 
convection to form inside the flux tube. 
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Table 5. Survey of numerical solutions obtained with dT/dz constant at the lower boundary and with parameters R = 10^, Co = 0.2, 
m = 1, 7 = 5/3, r = 3. The star superscript in the Q column indicates time dependent solutions. 



Q 


Q 


cr 


e 


Ro 


max|T| 


max(Mach) 


max \ U(j)\ 


max Br 


max 


max Bz 


range 





32 


1 


10 


oo 


1.2 


0.6 


— > 


6.3 


— > 


21.8 


(-7o,14zj 





128 


1 


10 


oo 


1.1 


0.4 


—> 


4.3 





13.4 


(-46,92) 





256 


1 


10 


oo 


1.1 


0.4 


— > 


3.7 


— > 


9.8 


(-38,72) 


0.1 


32 


0.1 


10 


23.7 


1.4 


1.3 


0.8 


6.3 


3.2 


16.9 


(-73,105) 


0.1 


32 


0.3 


10 


42.2 


1.2 


1.1 


0.5 


4.3 


3.0 


17.1 


(-56,106) 


0.1 


32 


0.6 


10 


59.8 


1.1 


0.6 


0.4 


5.0 


3.8 


31.0 


(-56,111) 


0.1 


32 


1 


10 


77.5 


1.1 


0.6 


0.5 


5.8 


3.7 


19.7 


(-59,149) 


0.1 


32 


1 


20 


205.5 


2.2 


0.8 


0.5 


4.3 


2.0 


14.3 


(-41,74) 


0.1 


128 


0.1 


10 


23.7 


1.4 


1.5 


0.5 


3.9 


1.5 


12.5 


(-38,80) 


n 1 
U. i 




U.o 


iU 


/I o o 


i.z 


U.o 


U.4 


Q Q 
O.O 


i.O 


1 1 o 
ii.z 


( Ad '7Q^ 

(-40, (o) 


0.1 


128 


0.6 


10 


59.8 


1.1 


0.6 


0.4 


4.1 


1.6 


11.4 


(-48,87) 


0.1 


128 


1 


10 


77.5 


1.1 


0.5 


0.4 


4.2 


1.6 


11.8 


(-4y,lUDj 


0.1 


128 


1 


20 


205.5 


1.9 


0.5 


0.3 


2.6 


0.8 


7.2 


(-25,40) 


n 1 

U. i 




n 1 


1 n 
iU 




1 Q 
i.O 


i. i 


n A 
U.4 


Q Q 
O.O 


1 n 
i.U 


O Q 

y.o 


(^-o4,dO ) 


0.1 


256 


0.3 


10 


42.2 


1.1 


0.6 


0.4 


3.2 


1.0 


8.3 


(-37,63) 


n 1 




U.o 


1 n 
iU 


oy.o 


i. i 


U.o 


n A 
U.4 


O.Z 


i. i 


Q Q 
O.O 


(^-o / ,oy ) 


0.1 


256 


1 


10 


77.5 


1.1 


0.4 


0.3 


3.3 


1.1 


8.4 


(-37,63) 


0.1 


256 


1 


20 


205.5 


1.6 


0.2 


0.1 


1.5 


0.4 


4.8 


(-15,29) 


0.2 


32 


1 


10 


38.7 


1.2 


0.5 


0.8 


3.3 


4.4 


21.5 


(-56,95) 


0.2 


32 


1 


20 


102.8 


2.2 


0.7 


0.7 


2.9 


3.0 


9.9 


(-44,54) 


0.2 


128 


1 


10 


38.7 


1.3 


0.4 


0.6 


2.8 


2.4 


12.2 


(-43,67) 


0.2 


128 


1 


20 


102.8 


1.9 


0.7 


0.5 


2.0 


1.4 


6.5 


(-24,35) 


0.2* 


256 


1 


10 


38.7 


(1.1,1.5) 


0.3 


0.5 


2.2 


(1.8,1.9) 


(5.6,12.8) 


(-35,57) 


0.2 


256 


1 


20 


102.8 


1.6 


0.2 


0.3 


1.4 


0.7 


4.3 


(-15,26) 


0.3 


32 


1 


10 


25.8 


1.2 


0.5 


0.9 


2.6 


4.8 


24.9 


(-54,110) 


0.3 


32 


1 


20 


68.5 


2.1 


0.6 


0.9 


2.2 


3.4 


10.2 


(-37,45) 


0.3* 


128 


1 


10 


25.8 


(1.1,1.3) 


(0.3,0.4) 


0.7 


(1.8,2.3) 


(2.6,2.8) 


(6.6,14.2) 


(-40,62) 


0.3 


128 


1 


20 


68.5 


1.9 


0.4 


0.6 


1.5 


1.8 


7.0 


(-24,27) 


0.3* 


256 


1 


10 


25.8 


1.1 


(0.2,0.3) 


0.5 


(1.4,1.8) 


(1.9,2.0) 


(6.8,7.4) 


(-33,46) 


0.3 


256 


1 


20 


68.5 


1.6 


0.2 


0.3 


1.1 


0.9 


4.1 


(-18,21) 



Increasing Q also brings time dependence to the solu- 
tion (Section I3.9p . For moderate Q values the weak con- 
vection cells oscillate horizontally inside the magnetic flux 
tube, while for large Q values the weak cells push period- 
ically through the edge of the flux tube into the field-free 
convection area. This time dependence can be reduced by 
increasing the strength of the magnetic field (Section l3.4|) . 

The collar flow around the magnetic flux tube is influ- 
enced by the strength of the convection and the temperature 
prescription at the lower boundary. By lowering the value of 
the Prandtl number (a), the convection becomes stronger 
and the size of the collar cell increases (Section 13. 5|) . The 
stronger convection pushes the magnetic flux tighter at the 
central axis so that the flux tube width decreases and the 
magnetic field strength on axis increases. For cr = 0.1 the 
collar flow survives a change of the lower boundary condition 
from a constant temperature to a constant dT/dz (Section 
13. op . However, for a ^ 0.3 the collar flow is destroyed and the 
magnetic fleld is dragged away from the central axis (Figure 
I17p . Weak convection cells form inside this wider flux tube. 

The azimuthal velocity and magnetic flelds are driven 
by the imposed Q, because in the absence of rotation these 
quantities have very small amplitudes, generated by the ini- 
tial plasma perturbation, which decay exponentially to zero 
with time (Section I3.ip . It follows that as Q increases, the 
magnitudes of and increase (Figure [TT]) . In contrast, 
the amplitudes of Ur and Uz hardly change with Q. For all 



values of Q the azimuthal flow pattern flts that of a Rankine 
vortex: in areas with strong magnetic fleld the plasma tends 
to rotate as a rigid body while around it a free vortex forms 
in the fleld-free convection areas. This means that max(i/^0) 
is located outside the flux tube edge. A flnite Q shortens 
the wavelength of convection in the radial direction, so that 
the initial convection annulus breaks up into more than one 
convection cell fSection l3.2p . The vortex forming around the 
flux tube is not dependent on the number of convection cells 
in the fleld-free region. 

The plasma inside the magnetic flux tube and the vortex 
around the tube flow prograde relative to the rotating cyUn- 
drical reference frame (Figure [S}. A retrograde or counter 
flow appears next to the outer wall of the cylinder. This 
counter flow is due to the fact that in our solution the ver- 
tical component of the angular momentum is zero relative 
to the rotating reference frame. We initialize the simula- 
tions with Lz = and the counter flow appears at the outer 
wall to maintain the sta tus quo. To obtain a retro grade flow 
at the central axis like[j ones Gallowavl ([l993»), we have 
to change the bottom boundary condition on the tempera- 
ture from constant T to constant dT/dz (Section I3.0p . This 
change in boundary condition also creates a strong horizon- 
tal magnetic component in the top layers of the numerical 
domain, which may rotate retrograde with the plasma at the 
axis and the outer wall (Figure I20p . Alternatively, by gen- 
erating weak turbulence inside the magnetic flux tube, it is 
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possible to perturb the rigid body rotation of the plasma 
inside the flux tube to such an extent that one gets pro- 
grade and retrograde flow inside the flux tube. This is more 
likely to happen with a constant dT/dz than for a constant 
T lower boundary condition. 

Unlike the azimuthal velocity, the azimuthal magnetic 
field is influenced by the structure of the convection cells. 
Max{B(f,) is confined to the strongest convection cell closest 
to the outer edge of the magnetic flux tube. For a constant 
T lower boundary condition this is usually the collar flow 
next to the magnetic flux tube. For constant dT/dz as lower 
boundary condition and a = 0.1, significant parts of form 
inside the weak convection in the flux tube as well as inside 
collar cell outside the tube (Figure [16]). For constant dT/dz 
and a ^ 0.3, the forms in the large convection cell around 
the flux tube, with a local maximum next to the flux tube 
f Figure fTTj) . The direction of B(p depends on the convection 
direction; for anticlockwise flow (as in the collar flow around 
the flux tube) it points in the negative direction and vice 
versa for clockwise flow. When the solution has one large 
convection cell with clockwise flow, which we obtain with a 
constant dT/dz at the bottom boundary, a local maxl^^l 
forms at the outer wall, but its radial width and amplitude 
is small so that it does not influence the numerical solution 
inside the domain (Section 13. 6|) . 

The current density in the (r, z) plane always forms 
around the local maxl^^l and flows in the same direction 
as the local convection. The azimuthal current density forms 
around the edge of the flux tube where the magnetic fleld 
lines have the largest gradients and curvature. This means 
any process that widens the flux tube, i.e. straightens the 
magnetic fleld lines, will decrease and vice versa. Increas- 
ing Q fSection [3.3p . the magnetic fleld strength fSection [3.4|) . 
the stratiflcation in the domain (Section 13. 7|) . and chang- 
ing the temperature lower boundary condition to constant 
dT/dz (Section 13. 6p lead to a decrease in the amplitude 
of while a lower Prandtl number (Figure [14]) increases 
the amplitude of j^. When the weak convection inside the 
magnetic flux tube becomes strong enough to bend the fleld 
lines, local maxima of Ij^I starts to form. 

Lowering the Prandtl number (a) increases the strength 
of convection (Section 13. 5p as well as thermal diffusivity. 
Thus stronger upflows lead to stronger localized heating in 
the upper layer, while the radial extent of the heated plasma 
is reduced (Figure [TJ]). The stronger convection also causes 
signiflcant variations in the density gradient inside the fleld- 
free convection area. In contrast, a flnite Q with a = 1 has 
little effect on the density inside the convection area (Section 
13. 3p . Only at the outer boundary does the rotation change 
the density gradient in a signiflcant way, but the radial ex- 
tent of this layer is small and does not influence the rest 
of the domain. Inside the magnetic flux tube the density is 
relatively unaffected for Q ^ 0.3. Relative large Q values are 
necessary to observe a signiflcant influence by the centrifugal 
force. 

To ascertain the effect of the lower boundary, we 
changed the temperature boundary condition (Section 13. 6p 
and the stratiflcation in the numerical domain f Section 13. 7|) . 
Increasing the stratification effectively increases the heat 
flux through the lower boundary into the domain. This 
widens the magnetic flux tube, allowing weak convection 
cells to form inside it. However, the convection in the fleld- 
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Figure 27. Velocity field on the (r,(f)) plane at z = 0.25 for 
Q = 32, (7 = l,r2 = 0.1, and constant T bottom boundary, shown 
in Figure [3] Arrows represent Ur and and colour the sound 
speed perturbation Cg. Max(cs)=0.54 is the light shade between 
the two convection cells at r = 1.9 and min(cs)=-0.03 the dark 
shade on the edge of the magnetic fiux tube at r = 0.5. 



free regions and the configuration of the magnetic field stay 
essentially the same (Table [3] and Figures [23] and [24]) . In 
contrast, changing the temperature prescription from a con- 
stant temperature to constant dT/dz drastically affected the 
solution. The bottom temperature reduces slightly (Figure 
[22]) . but this does not account for the changes in the solu- 
tion. The amplitude of the convection reduces significantly 
(Table [2} and for a ^ 0.3 the fiow pattern and magnetic 
field configuration change radically (Figures [T6l and [T7]) . 

The numerical solutions obtained in this study point to 
a specific radial profile for azimuthal velocities in sunspots 
that rotate around their own axis. Inside the umbra, where 
the vertical magnetic field component is strong, the plasma 
rotates as a rigid body while the convection around the um- 
bra is in the form of a vortex. This profile is supported by 
most of the observations. Photospheric observations place 
the maximum azimuthal velocity inside the penumbra, while 
helioseismic observations show a vortex fiowing around the 
fiux tube in the convection zone. The typical azimuthal ve- 
locity (us) in the p h otosphere is of the o r der 10 ~^ km s~^ 
(|Brown et al ] |2003l ). IZhao k KosovichevI (|2003l ) measured 
max(i^0) ^ 0.5 km s~^ below the photosphere at depths 
0-3 and 9-12 Mm. This compares well with our measured 
velocities in Tables [J] and [5] where for low angular velocity 
{VL — 0.1) we obtain max(n0) ^ O(10~^) km s~"^, taking a 
sound speed of 1.29 km s~^ as reference speed. 

We present Figures [27l and [28l to facilitate comparison of 
our results with l ocal helioseismic m easurements, of which 
Figures 6 to 8 in iKosovichevI (|2002l ) are examples. Figure 
[271 shows the fiow in the (r, (j)) plane for a constant T and 
Figure [28] for a constant dT/dz bottom boundary. These 
planes correspond to depth z = 0.25 in Figures [3l and [TTl re- 
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Figure 28. Velocity field on the (r, ^) plane at z = 0.25 for 
Q = 32, a = 1, Q = 0.1, and constant dT/dz bottom boundary, 
shown in Figure fTTl The diagnostics has the same meaning as in 
Figure [27] Max(cs)=0.04 is the light shade around the magnetic 
flux tube at r = 1 and min(cs)=-0.12 the dark shade at the outer 
wall at r = 3. 



spectively, so that Figure [27] represents two convection cells 
with a collar flow around a well-deflned magnetic flux tube, 
while Figure [28] represents one outflowing convection cell 
that drags the magnetic fleld lines away from the central 
axis. The arrows show that Ur dominates for Q = 0.1, 
with azimuthal flow patterns more visible in the inner radius 
closer to the magnetic flux tube. The size of Ucf, relative to Ur 
increases with Q (Figure [TT]) . so that the Rankine vortex be- 
comes more visible for higher values of ^1. In Figures [27] and 
[28] the outer boundary condition Ur = still holds, with ar- 
rows chosen close to the boundary showing that the flow has 
finite size next to the boundary. The colour palette shows 
the perturbed sound speed (cs)in the plane. Where there is 
an upwelling the plasma is heated and vice versa. Figure [27] 
shows the warmer plasma - and hence higher sound speed - 
between the two convection cells and Figure [28] at the upflow 
next to the magnetic flux tube. Downflow with its accompa- 
nied cooler plasma - and hence lower sound speed - occurs 
around the edge of the flux tube for Figure [27] and at the 
outer wall for Figure 1281 The flux tube itself is cooler than 
the rest of the surrounding plasma (Figure [27]), but where 
weak convection inside it exists, it starts to heat up (Figure 
I28p . The difference in max (cs) between Figures [27] and [28] is 
thus due to the radical different flow pattern for each case. 

We generate azimuthal flow in this axisymmetric model 
by rotating the cyUnder around its axis at a constant angular 
velocity. As a result we obtain time independent solutions, 
in contrast to the highly time dependent observations. For 
low angular velocities the flow inside the magnetic flux tube 
and the vortex flow are prograde. Due to the fact that our 
model conserves the vertical component of the angular mo- 
mentum, a retrograde flow appears next to the outer wall. 



We find that high angular velocities tend to break the um- 
bra into concentric rings and introduce time dependence in 
the form of periodic behaviour in the radial direction. These 
phenomena have not been observed in sunspots and may be 
due to the axisymmetry in our model. It is more likely that 
our numerical results obtained with low angular velocities 
are reaUstic models of solar observations. 
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